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Abstract 

We implement a high-order numerical scheme for the entropy-based moment closure, the so-called Mjv model, 
for linear kinetic equations in slab geometry. A discontinuous Galerkin (DG) scheme in space along with a 
strong-stability preserving Runge-Kutta time integrator is a natural choice to achieve a third-order scheme, 
but so far, the challenge for such a scheme in this context is the implementation of a linear scaling limiter 
when the numerical solution leaves the set of realizable moments (that is, those moments associated with a 
positive underlying distribution). The difficulty for such a limiter lies in the computation of the intersection 
of a ray with the set of realizable moments. We avoid this computation by using quadrature to generate 
a convex polytope which approximates this set. The half-space representation of this polytope is used to 
compute an approximation of the required intersection straightforwardly, and with this limiter in hand, the 
rest of the DG scheme is constructed using standard techniques. We consider the resulting numerical scheme 
on a new manufactured solution and standard benchmark problems for both traditional Mjv models and the 
so-called mixed-moment models. The manufactured solution allows us to observe the expected convergence 
rates and explore the effects of the regularization in the optimization. 
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1. Introduction 

Moment closures are a class of spectral methods used in the context of kinetic transport equations. An 
infinite set of moment equations is defined by taking velocity- or phase-space averages with respect to some 
basis of the velocity space. A reduced description of the kinetic density is then achieved by truncating this 
hierarchy of equations at some finite order. The remaining equations however inevitably require information 
from the equations which were removed. The specification of this information, the so-called moment closure 
problem, distinguishes different moment methods. In the context of linear radiative transport, the standard 
spectral method is commonly referred to as the P n closure [51] , where N is the order of the highest-order 
moments in the model. The Pat method is powerful and simple to implement, but does not take into account 
the fact that the original function to be approximated, the kinetic density, must be non-negative. Thus Pat 
solutions can contain negative values for the local densities of particles, rendering the solution physically 
meaningless. 

Entropy-based moment closures, referred to as Mat models in the context of radiative transport Enm, have 
all the properties one would desire in a moment method, namely positivity of the underlying kinetic density^ 


^ Positivity is actually not gained for every entropy-based moment closure but is indeed a property of those models derived 
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hyperbolicity of the closed system of equations, and entropy dissipation m- Practical implementation of 
these models has been traditionally considered too expensive because they require the numerical solution of 
an optimization problem at every point on the space-time grid, but recently there has been renewed interest 
in the models due to their inherent parallelizability |12j . However, while their parallelizability goes a long 
way in making M^v models computationally competitive, in order to make these methods truly competitive 
with more basic discretizations, the gains in efficiency that come from higher-order methods will likely be 
necessary. Here the issue of realizability becomes a stumbling block. 

The property of positivity implies that the system of moment equations only evolves on the set of so-called 
realizable moments. Realizable moments are simply those moments associated with positive densities, and 
the set of these moments forms a convex cone which is a strict subset of all moment vectors. This property, 
while indeed desirable since it is consistent with the original kinetic density, can cause problems for numerical 
methods. Standard high-order numerical solutions to the Euler equations, which indeed are an entropy-based 
moment closure, have been observed to have negative local densities and pressures [35]. This is exactly loss 
of realizability. 

A recently popular high-order method for hyperbolic systems is the Runge-Kutta discontinuous Galerkin 
(RKDG) method |5l|6j. An RKDG method for moment closures can handle the loss of realizability through 
the use of a realizability (or “positivity-preserving”) limiter [35] , but so far these have been implemented 
for low-order moment systems (that is A = 1 or 2) [23| because here one can rely on the simplicity of the 
structure of the realizable set for low-order moments. For higher-order moments, the realizable set has complex 
nonlinear boundaries: when the velocity domain is one-dimensional, the realizable set is characterized by the 
positive-definiteness of Hankel matrices IH1I2H1; in higher dimensions, the realizable set is not well-understood. 
In [1], however, the authors noticed that a quadrature-based approximation of the realizable set is a convex 
polytope. With this simpler form, one can now actually generalize the realizability limiters of ^51155] for 
moment systems of (in principle) arbitrary moment order. Furthermore, this approximation of the realizable 
set holds in any dimension. 

In this work we begin by reviewing our kinetic equation, its entropy-based moment closure, and the concept 
of realizability in Section Then in Section [^ we outline how we apply the Runge-Kutta discontinuous 
Galerkin scheme to the moment equations. Here the key ingredients are a strong-stability preserving 
Runge-Kutta method, a numerical optimization algorithm to compute the flux terms, a slope limiter, a 
realizability-preserving property for the cell means, and the realizability limiter. In Section]^ we present 
numerical results using a manufactured solution to perform a convergence test, as well as simulations of 
standard benchmark problems. Finally in Section]^ we draw conclusions and suggest directions for future 
work. 


2. A linear kinetic equation and moment closures 

We begin with the linear kinetic equation we will use to test our algorithm and a brief introduction to 
entropy-based moment closures and the concept of realizability. More background can be found for example 
in miiiniiiii and references therein. 

2.1. A linear kinetic equation 

We consider the following one-dimensional linear kinetic equation for the kinetic density '0 = tjj(t,x,fj,) > 0 
in slab geometry, for time t > 0, spatial coordinate x G X = (xljXr) C M, and angle variable /i G [—1,1]: 

dt'tp + = asC {^) + S, (2.1) 


from important, physically relevant entropies. 
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where era are CTs are the absorption and scattering interaction coefheients, respectively, which throughout 
the paper we assume for simplicity to be constantsj^and S a source. The operator C is a collision operator, 
which in this paper we assume to be linear and have the form 

1 .1 

dfi'. ( 2 . 2 ) 

1 J-i 

We assume that the kernel T is strictly positive and normalized to , fj,)dii' = 1. A typical example 

is isotropic scattering, where = 1 / 2 . 

Equation (|2.1[) is supplemented by initial and boundary conditions: 


p{t,XL,fJ') = 

= V’L(i,/i) , 

t > 0 , 

/r > 0 , 

(2.3a) 

p{t,x^,p) = 


t > 0 , 

77 < 0 , 

(2.3b) 

p{t),x,p) = 

= 'pt=o(,X,p) , 

X G (xL,a:R), 

p G [—1,1] , 

(2.3c) 


where V'Lj V'Rj and ipt=o are given. 



2.2. Moment equations and entropy-based closures 


Moments are defined by angular averages against a set of basis functions. We use the following notation for 
angular integrals: 


(</)=/ 1 


for any integrable function (p = (/(/r); and therefore if we collect the basis functions into a vector b = b(/r) = 
{bo{n), bi{fi ),..., bffiti))'^, then the moments of a kinetic density (p = p{p,) are given by u = (b^). 

A system of partial differential equations for moments u = u(t, x) approximating the moments {h^p) (for 
the Ip which satisfies ( 2 . 1 )) can be obtained by multiplying ( 2 . 1 ) by b, integrating over p, and closing the 
resulting system of equations by replacing ip where necessary with an ansatz ^p^^ which satisfies u = {h^p^^) . 
The resulting system has the form [HHO] 


where 


dtu + a,,f(u) + (TaU = crsr(u) + (hS ), 


f(u) := [fih'ipu) and r(u) := (bC 


(2.4) 


In an entropy-based closure (commonly referred in standard polynomial bases as the M^v model), the ansatz 
is the solution to the constrained optimization problem 


Tpu = argmin{(r 7 ((/)) : (hp) = u} , 


(2.5) 


where the kinetic entropy density r] is strictly convex and the minimum is simply taken over functions 


p = p{ti) such that {rjiP)) is well defined. The optimization problem (2.5) is typically numerically solved 
through its strictly convex finite-dimensional dual. 


i:(u) := argmin (77*(b^Q;)) — 




( 2 . 6 ) 


where 77 * is the Legendre dual of rj. The first-order necessary conditions for the dual problem show that the 

(2.7) 


solution to the primal problem (2.5) has the form 

■pu = vl (b'^Q:(u)) 


^ All results here can be generalized to spatially inhomogeneous interaction coefficients. 
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where ry' is the derivative of 77*. 


The entropy rj can be chosen according to the physics being modeled. As in [12] we use Maxwell-Boltzmann 
entropj0 


r]{z) = z log(z) — z. 


For the Maxwell-Boltzmann entropy (2.8), 77* (y) = 77^(2/) = exp(7/). 

In this paper we consider both the monomial moments, dehned by the basis 

p(m) := 


( 2 . 8 ) 


and the so-called mixed-moments m which contain the usual zeroth-order moment but half moments in 
the higher orders. This is achieved using the basis functions 


m(/r) ;= (l, 77+, 77 .,..., 77 +, 77^)^ , 

where 77+ = max(77,0) and fi- = min(77, o)Q Mixed-moment models, which we refer to as MM^r, have been 
introduced to address disadvantages like the zero net-flux-problem and unphysical shocks in full-moment 

models [iniiizj- 

We close this section by quickly noting that the classical Pn approximation m is an entropy-based moment 
closure by choosing the basis b as the Legendre polynomials and using the entropy 

77(z) = 


This results in the ansatz 


i/ju = b^ a. 


which clearly is not necessarily positive. Nonetheless the resulting moment system is linear, simple to 


compute, and for high values of N provides good baseline solutions to the original kinetic equation (2.1). 


2.3. Moment realizability 


Since the underlying kinetic density we are trying to approximate is nonnegative, a moment vector only 
makes sense physically if it can be associated with a nonnegative density. In this case the moment vector is 


called realizable. Additionally, since the entropy ansatz has the form (2.7), in the Maxwell-Boltzmann case 


the optimization problem (2.5) only has a solution if the moment vector lies in the ansatz space 


A := {(bexp (b^a)) : a G . 


In our case, where the domain of angular integration is bounded, the ansatz space A is exactly equal to the 
set of realizable moment vectors M- Therefore we can focus simply on realizable moments, so in this section 
we quickly review their characterization in the cases of exact and approximate integration. 


^ Indeed in a linear setting such as ours, any convex entropy 77 is dissipated by ( |2.1[ l, so we have some freedom. We focus on 
the Maxwell-Boltzmann entropy because it is physically relevant for many problems, gives a positive ansatz and also allows 
us to explore some of the challenges of numerically simulating entropy-based moment closures. 

Notice that in the mixed-moment case, there are 2^” -h 1 basis functions instead of N" -|- 1 as in the monomial case. However, 
for clarity of exposition, for most of the paper we will assume b has -h 1 components, though everything applies to the 
mixed-moment case as well. 
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2.3.1. Classical theory 


Definition 2.1. The realizable set TZ-^ is 

= {u : 30(/r) > 0, (</)) > 0, such that u = (b^)} . 

Any (j) such that u = (bi^) is called a representing density. 

The realizable set is a convex cone. 

In the monomial basis b = p, a moment vector is realizable if and only if its corresponding Hankel matrices 
are positive definite [28j . When a moment vector sits exactly on dlZ^, there is only one representing density, 
and it is a linear combination of point masses [S]. In this case, the corresponding Hankel matrices are 
singular. This also causes the optimization problem to be arbitrarily poorly conditioned as the moment 
vector approaches dTZp [5]. 

Realizability conditions in the mixed-moment basis b = m are given in |27j again using Hankel matrices for 
each half-interval [—1,0] and [0,1] as well as another condition to “glue” the half-interval conditions together. 
In this case, only a subset of the moment vectors on dTZ^ have unique representing densities, but those that 
do include point masses. 


2.3.2. The numerically realizable set 

In general, angular integrals cannot be computed analytically. We define a quadrature for functions 
(j): [—1,1] —>■ K by nodes and weights such that 

UQ 

« {(j)) 


Below we often abuse notation and write {(j)) when in implementation we mean its approximation by 
quadrature. Then, as defined in [T], the numerically realizable set is 


UQ 


= < u : 3/j > 0 s.t. u = y] Wih{iJ.i)fi 


i=l 


Indeed, when replacing the integrals in the optimization problem (2.5) with quadrature, a minimizer can 
only exist when u S 7^®. It is straightforward to show that, as expected, 7^® C TZ^. 


The numerically realizable set 77.® is the convex cone generated by L -i’ normalized moment 

vectors: 

'^b\uo=i •= G 77^ : Mo = 1} ; 

and Lq-i interior of a convex polytope: 

Proposition 2.2 ([I]). For any quadrature Q with positive weights Wi, and for simplicity assuming 6o(m) = 1? 

=intco{b(^i)}r®i 

where int indicates the interior and co indicates the convex hull. 
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3. Realizability-preserving discontinuous Galerkin scheme 


In this section we introduce our high-order numerical method to simulate the moment system (2.4). We use 


the Runge-Kutta discontinuous Galerkin (RKDG) approach [5l [6] and recent techniques for the numerical 
solution of the defining optimization problem [1] (2.5)-( [2^ . Finally in this section we discuss the crucial 
issue of realizability and our linear scaling limiter to handle non-realizable moments in the solution. 


3.1. The discontinuous Galerkin formulation 

We briefly recall the discontinuous Galerkin method for a general system with source term: 

5tu-ha2;f(u) = s(u), 


(3.1) 


where in our case s(u) := tTsr(u) — -|- (bS'). We follow the approach outlined in a series of papers by 

Cockburn and Shu We divide the spatial domain (xljXr) into J cells Ij = {xj_i/2,Xj+ii2), where the 

cell edges are given by Xj^xji = Xj ± Aa;/2 for cell centers Xj = xl + {j — l/2)Ax, and Ax = (xr — x\)lJ . 
For each t, we seek approximate solutions Uh{t,x) in the finite element space 


Vl: = {v€ L1(xl, Xr) : v\i^ e P'=(/,) for j S {1,..., J}}. 


(3.2) 


where P^{I) is the set of polynomials of degree at most k on the interval I. We follow the Galerkin approach: 


replace u in (3.1) by a solution of the form U/j £ then multiply the resulting equation by basis functions 


Vh of and integrate over cell Ij to obtain: 

dt / U/i(t,x)nft,(x)dx-hf(u/,(t,x“^^/2))'^'*(^7+1/2) “ x^_^^.^))vhix+_j^^.^) 

- / i{\ih(t,x))d^Vh{x)dx = / s{\ih{t,x))vh{x)dx (3.3a) 

dij dij 

/ UhiO,x)vh{x)dx = / Ut=o{x)vh{x) dx (3.3b) 

Ji, Ji, 

where ^j±i /2 denote the limits from left and right, respectively, and Ut^o is the initial condition. 

In order to approximately solve the Riemann problem at the cell-interfaces, the fluxes f(u?i(t, x^^j^^ 2 )) tii® 
points of discontinuity are both replaced by a numerical flux f (u;i(t, x“^j^y 2 )i 2 ;^+i/ 2 )); thus coupling 
the elements with their neighbors [3T]. In this paper we use the global Lax-Friedrichs flux0 

f(v, w) = i (f(v) -t- f(w) - C(w - v)), 

The numerical viscosity constant C is taken as the global estimate of the absolute value of the largest 
eigenvalue of the Jacobian df/du. We use (7=1, because for the moment systems used here it is known 
that the largest eigenvalue is bounded by one in absolute value: 

Lemma 3.1. The eigenvalues of the Jacobian di/du are bounded in absolute value by one. 


® The local Lax-Friedrichs flux could be used instead. This would require computing the eigenvalues of the Jacobian in every 
space-time cell to adjust the value of the numerical viscosity constant C but would possibly decrease the overall diffusivity of 
the scheme. However, since we are considering high-order methods, the decrease in diffusivity achieved by switching to the local 
Lax-Fridrich flux should be negligible. 
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Proof. For convenience, we present a slight generalization of Lemma 4 in |23j . 

We define J(q!) := (/rbb^? 7 "(b^a)) and H(a) := (bb^ 7 y"(b^Q:)). Using the properties of H(a) given in 
[20) . by applying the chain rule we have 

^ = J(a(u))^||^ = J(«(u))H(«(u))-i. (3.4) 

If J(Q;(u))H(a(u))“^c = Ac for some c 7 ^ 0, then for d = H(q:(u))“^c we also have J(Q:(u))d = AH(Q:(u))d, 
and thus 

^ d^J(d(u))d _ (M(b^d)%:(d(u))) ^ ^ 
d^H(Q:(u))d (^(b^d)^ ry"(d(u))^ 

The last inequality follows from the facts that |/i| < 1 and that 77 " > 0, the latter of which is a consequence 
of the strict convexity of 77 . □ 


On each interval, the DG approximate solution u/j can be written as 

k 

i=0 


X — Xj 

Ax 


(3.5) 


where {tpo, Pi, ■ ■ ■, Pk} denote a basis for P^{[—l/2, 1/2]). It is convenient to choose an orthogonal basis, so 
we use Legendre polynomials scaled to the interval [—1/2,1/2]: 

po{y) = l, (fi{y) = 2y, ip 2 iy) = ^{I2y'^ - 1), .. . 

With an orthogonal basis the cell means u^- are easily available from the expansion coefficients Uji 


Ic 

1 /* ^ ^ 

■= aU I. 


Ax 


p^ 


i=0 


X — Xj 

Ax 


dx = u°(t) 


We collect the coefficients iP (t) into the (A: + 1) x (N + 1) matrix 


Ui(t) = 




Using the form of the approximate solution in (3.5), we can write (^3.3b in matrix form: 


MdtUj + F(uj_i,U 7 ,Uj+i) - V(uj) = S(uj) 

(Muj(0))^^= / ue^t=o{x)pi 

Ji, 

for / e {1, 2,..., J}, with 


{M)u = Pi 


Ax 


Pi 


Ax 


Ax 


dx. 


dx 


(V(uj(t)))if fe{uj{t,x))d,,<fi dx, 

{S{Uj{t)))u = Si,{Uj{t,x))(fi 


Ax 


dx, 


(3.6a) 

(3.6b) 

(3.7a) 

(3.7b) 

(3.7c) 

(3.7d) 

(3.7e) 
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where U£^t=o, fi, fe, and se are the ^-th components of Ut= 0 j fi f; and s respectively. Notice that M is 
diagonalized by the choice of an orthogonal basis {‘fi}- We can write (3.6aI as the ordinary differential 
equation 


dtxij = LhiUj-i,Uj,Uj+i), for j G {1,..., J} and t S (0,r), 


(3.8) 


with initial condition specified in (3.6b). 


Boundary conditions are incorporated into the quantities Uo(t, 0 : 1 / 2 ) and a:j+i/ 2 ), which we have not 

defined yet but appear in the numerical flux (3.7b) for the first and last cells. To define these terms, first we 


smoothly extend the definitions of and fi) to all fi (note that while moments are defined using 

integrals over all fi, the boundary conditions are in (2.3a)-( '2.3b[ ) only defined for fi corresponding to incoming 
data)nthen simply take Uo(t, 0 : 1 / 2 ) := {hipi,{t,-)) and uj+i(t,o:j+i/ 2 ) := (b^R(t, •)). This completes the 
spatial discretization. 


In this paper, except for some of the convergence tests in Section 4.1.1 we use quadratic polynomials (k = 2) 
resulting in a third-order approximation. The integrals in (3.7) are computed using quadrature exact for 


polynomials of degree five to ensure the numerical scheme is third-order convergent. We use the four-point 
Gauss-Lobatto quadrature rule since the function evaluations at the interval boundaries can be reused for 
the numerical fluxes F. 


3.2. Runge-Kutta time integration 


For a fully third-order method, we require a time-stepping scheme for (3.8) that is at least third-order 


We use the standard explicit SSP(3, 3) third-order strong stability-preserving (SSP) Runge-Kutta time 
discretization introduced in [5S]. Let {t^y^Lo denote time instants in [0,tt] with = nAt, and for each cell 
j G {1,2,..., J} let the initial coefficients u/(0) be defined as in (3.6b). Then for n G (0,1,..., W — 1} we 
compute u/(t”) as follows: 

= u/(t”) -h AtL/j(u/_i(t”),U/(t”),Uj+it”)); 


u 


( 2 ) 


This specific Runge-Kutta method is a convex combination of forward Euler steps, a property which below 
helps us prove that the cell means of the internal stages are realizable. 

A scheme of higher order could be achieved by increasing the degree k of the approximation space as well 
as the order of the Runge-Kutta integrator. Unfortunately SSP-RK schemes with positive weights can at 
most be fourth order mm- A popular solution is given by the so-called Hybrid Multistep-Runge-Kutta 
SSP methods. A famous method is the seventh-order hybrid method in m while recently two-step and 
general multi-step SSP-methods of high order have been investigated mm- 


(1) ,■',(1) -(i) 


3 . 3 . Numerical optimization 

In order to evaluate f(u) and r(u) on the spatial quadrature points in each Runge-Kutta stage, we first 


compute the multipliers Q:(u) solving the dual problem (2.6). For the Maxwell-Boltzmann entropy the dual 
objective function and its gradient are 


f{a) := (exp(b^a)) — and g(Q:) = (bexp(b^a)) — u. 


® Although this is indeed the most commonly used approach, its inconsistency with the original boundary conditions 
(|2.3a[|-|2.3b[l, is still an open research topic 11714191 [24l 1301 . 
























respectively. 

We use the numerical optimization techniques proposed in [1] . The stopping criterion for the optimizer is 
given by 

I|g(a)ll2 < T, 

where || • II 2 is the Euclidean norm, and t is a user-specified tolerance, and we also use the isotropic 
regularization technique to return multipliers for nearby moments when the optimizer failsj^ Isotropically 
regularized moments are defined by the convex combination 

v(u, r) := (1 - r)u + ruoUiso, 

where 

Uiso = 2 (^) (3-9) 

is the moment vector of the isotropic density </>(/i) = 1/2. The form of v(u,r) is also chosen so that v(u, r) 
has the same zeroth-order moment as u. We define for an outer loop an increasing sequence {r^} for 
m = 0,1,..., We begin at m = 0 with r = tq ;= 0 and only increment m if the optimizer fails to 

converge for v(u, r^) after iterations. It is assumed that is chosen large enough that the optimizer 

will always converge for v(u,rm„,,x) for any realizable u. 


3.4- Realizability preservation and limiting 


In order to evaluate the flux-term f(u_,(t, Xjq)) at the spatial quadrature nodes Xjq in the /-th cell, we at least 
need Uj(t,Xjq) € TZ^ for each node, although when the angular integrals are approximated by quadrature, 
we in fact need Uj(t,Xjq) G 7Z^. Unfortunately higher-order schemes typically cannot guarantee this, as 
has been observed in the context of the compressible Euler equations (which are indeed in the hierarchy of 
minimum-entropy models) in [115] . 


We can, however, Hrst show that, when the moments at the quadrature nodes are realizable, our DG scheme 
preserves realizability of the cell means Uj(t) under a CFL-type condition. With realizable cell means 
available, we then apply a linear scaling limiter to each cell pushing Uj(t,Xjq) towards the cell mean and 
thus into the realizable set for each node Xjq. 


Following the arguments in [33l [3l] , this limiter does not destroy the accuracy of the scheme in case of 
smooth solutions if u,- is not on the boundary of the realizable set. We test this numerically in Section 4.1.2 


3.4-1- Realizability preservation of the cell means 

To prove realizability preservation of the cell means we will need three main ingredients: first, an exact 
quadrature to represent the cell means using point values from the cell; second, a representation of the 
moments collision operator; and finally a lemma that allows us to add the flux term without leaving the 
realizable set. 

First, following [3511371, we consider the Q-point Gauss-Lobatto rule and use its exactness for polynomials of 
degree k < 2Q — 3 to write the cell means as 




If ^ 

— Uj{t,x)dx = '^WqUj{t,Xjq), 

“'b q=l 


(3.10) 


^ The optimizer can fail for two reasons: either the Cholesky factorization required to find the Newton direction fails or the 
number of iterations reaches a user-specified maximum kmax- 
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where Xjq € Ij are the quadrature nodes on the j-th cell, and Wq are the weights for the quadrature on the 
reference cell [—1/2,1/2]. We choose the Gauss-Lobatto quadrature in particular because it includes the 
endpoints, that is 

^7+1/2 = XjQ (3.11) 


xt 1 /r, = x^i and a;,., 


3 - 1/2 


Secondly, we note that since the collision kernel T in (2.21 is positive by assumption, the moments of the 
collision operator applied to the entropy ansatz can be written as the difference between a realizable moment 
vector uc(u) and the given moment vector u: 


r(u) = (bC 


= uc(u) - u. 


(3.12) 


For example, in the case of isotropic scattering (T = 1/2), we have uc(u) = UoUjso (where Ujso are the 
moments of the normalized isotropic density, see (3.9)). 

Finally, we use the following lemma: 

Lemma 3.2. // u G and a — |5| > 0, then au + 6f(u) G TZ^^. 


Proof. Since u G TZ^, then there exists a solution to the minimum entropy problem so that u = (b^i/u) 
and f(u) = {fihtpu). Thus 

au + 5f(u) = ^b (a + bfj.) 'i/u^ • 

Since ^ G [—1,1], the affine polynomial satisfies a + bfj. > a — \b\, so hy assumption we have (a + &^)V’u > 0, 
which is a nonnegative measure representing au + bf (u). □ 

Theorem 3.3. Assume S > 0 and 2Q — 3 > fc, and let at := as + a^- Consider the cell means at time 
instants C, 




If the moment vectors := Uj{C,Xjq) at each spatial quadrature point Xjq are in TZ^^ (or TZ^) and the 
CFL condition 


At 

— < n;Q(l - CTtAt). 


(3.13) 


holds, where wq denotes the quadrature weight of the last quadrature weight of the Q-point Gauss-Lobatto 

hei 

time step for (3.6) are also in IZ^ (respeetively TZ^). 


quadrature rule on the referenee eell [—1/2,1/2], then the eell means computed by taking a forward-Euler 


Proof. The following arguments only use the fact that the realizable set is a convex cone, and therefore can 
be applied to either TZ^, or 7^® in exactly the same way. For clarity of exposition, we also begin with the 
case S' = 0. 


By the orthogonality of our basis, Uj{t) = u^^\t). Therefore we use the subequations in (3.6a) for u 
(where, in particular we have = 0). We use the notation Uj{x) = Uj(t", x), so that with a forward-Euler 
approximation for the time derivative (3.6a) gives 


( 0 ) 


-n+1 _ 


= u/ + At - 




Ax 


- CTaU” crsr(u/) . 
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Now using the fact that the cell interfaces are the first and last quadrature nodes (recall (3.11)), we now 


substitute the definition of f and the appropriate representation of the moments of the collision operator 


(3.12), then use the quadrature formula for the cell means (3.10), and finally collect terms to get: 

1 


u”+i =u” + At ( - 


2Ai (' ' (“"j+i)i) ■ (“5+1)1 ■ 

(f(“5-i)Q)+f(“5i)-(“?i-u;‘_„o))) 


- (TaU” + (Ts(uc(u") - U") 

= E (f + f ("0+i)i) - ("0+1)1 - 

- (f ("O-Dq) + f «i) - ("fi - "0-i)Q 
Q 

-CTt'^ WqUjg + tTgUc (u 


9=1 


Q-1 


= ^ Wq{l - Atat)u^q + AtusUc (u”) 

q=2 

+ ("0+1)1 - f ("0+i)i) + "0-i)Q + f ("0 -i)q)) 

At 


Wi 




(3.14a) 

(3.14b) 

(3.14c) 

(3.14d) 


Keeping in mind that we have assumed that each moment vector is realizable, we consider each of the 
final lines: 


• If (TtAt < 1, a condition which is indeed weaker than (3.13), the expression in the first line, (3.14al 


is a positive linear combination of realizable moments. Since the realizable set is a convex cone, this 
expression is realizable. 


• The terms in the second line (3.14b) can be shown to be realizable by two applications of Lemma 3.2 
with 0 = 1 and b = ±1. 


(3.13) (and recalling that wi = wq). 


• The expressions in the last line, (3.14cl and (3.14dl, are each realizable according to Lemma 3.2 and 


Finally, 0"+^ is realizable since it is a sum of realizable moment vectors. 


When S' > 0, notice that this simply adds to (3.14) the term (bS), which is realizable and thus does not 
affect the conclusion. □ 


Since we are using SSP-Runge-Kutta time-stepping schemes, whose stages are convex combinations of forward 
Euler steps. Theorem |3.3| guarantees that under the appropriate CFL condition, the cell means for every 
Runge-Kutta stage are realizable. In particular, the SSP(3, 3) which we use is a convex combination of Euler 
steps all with time step At. 
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(a) Before limiting; notice the nonrealizable moments at 
the end points. 



Figure 1: Application of the realizability preserving limiter to the Mi system for a quadratic polynomial. Here, (uo,ni) are 
realizable only if |ui| < no. Black squares indicate the spatial quadrature points Xq. 


S.^.2. Realizability-preserving limiter 

Theorem |3.3| makes the assumption that the point-values of the local DG polynomials are realizable, and this 
must also hold in order to evaluate the flux f at the moment vectors on the quadrature points. This can be 
achieved by applying a linear scaling limiter in each cell. Recall the deflnition of Uj(t,x): 

i=0 ^ ' i=l ^ ' 

We can see that using convexity of the realizable set, if u^- is realizable, then for each quadrature point there 

exists a 0 e [0,1] such that 

:= OUjit) + (1 - 0)Uj(t,Xjq) 

is realizable. Indeed, by inserting the dehnition of Uj(t,Xjq) from above, we can write the limited moment 
vector as 


u“ 


{t,Xjq) = Uq{t) + (1 - 0) 


k 

E-f 

2 = 1 




^31 


Ax 


thus when limiting is necessary, the higher-order coefficients are damped while the cell mean remains 
unchanged. An example of the limiting process is illustrated in Figure]^ which considers the following 
polynomial representation of an Mi solution: 



f 1 0.5 

1 ^ 0.8 0.2 





After limiting with 6 > 9/11 the vector (mq,Mi) becomes realizable. 
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Since the full realizable set TZ^ is characterized by the positive-definiteness of the associated Hankel matrices, 
computing the smallest 0 such that Uj{t,Xjq) G is in general difhcult. Furthermore, for higher-dimensional 
problems (when dimension of the angular domain is more than one) the realizable set is in general not 
well-understood. 

This computation is however easier for the numerically realizable set 72,® using its half-space representation. 
This representation is also intriguing because it extends to higher-dimensional problems, since even in these 
cases 72.® remains the interior of a cone generated by a convex polytope. 

In the following for clarity of exposition we omit the time arguments and spatial-cell indices, therefore using 
u to indicate the (always realizable) moment vector at the cell mean and Ug to indicate the (not necessarily 
realizable) moment vector at a quadrature point. In an implementation, the limiter is applied to every 
quadrature point in every cell. 

To discuss the computation of 6 , first we assume that the moment vectors u and have been scaled such 
that max(uo, Wqo) = 1/2, where uq and Uqo are the zero-th components of u and Ug, respectively!^ Then the 
limited moment vector u® := du -|- (1 — 0)ug satisfies < I for all 6 G [0,1], and without loss of generality 
we can apply the limiter to move the moments into the bounded set 

^bL„<i = {ue^b |uo<i}. 


instead of the full, unbounded set 7^®. Now, using that 7^® is the cone generate by "7^® L^-i Proposition 
can be written as the interior of a convex polytope 


2.2 




«0<1 


7?® 

'^'b 


uo<l 


= |Au : A e (0,1) and u S "7^® Lo=i} = int co {0, b(/ri), b(/r 2 ), ■ •. ,b(^„g)} . 


As a convex polytope, it has a half-space representation [3] of the form 

^b Lg<i = {u e : afu < h,iG {!,...,d}} , 

where d is the number of facets of the polytope. For each i-th facet, the vectors a.^ and scalars are computed 
from the set of vertices defining the facet. These sets of vertices can be computed using standard convex-hull 
algorithms, and for our implementation we used the Matlab routine convhulln. These are fixed throughout 
the computation and therefore can be precomputed. 

Candidate values 9qi for the 9 which ensures u® G can now be computed for each i-th facet by 


a^ F (1 dg2)llg) - 


_ b,-afuq 

^ni — 


af (u - Ug) 


If we wanted to limit exactly to 972®, we would first take 


0 

max{6»gi : 9qi G [0,1]} 


if there is no 9qi G [0,1], 
else; 


and then for the cell in question, we would apply the largest 9q from the quadrature nodes: 9^^ := max{dg^}. 
This would ensure that the moment vectors at each quadrature node in the cell are in the realizable set or on 
its boundary. 

In practice, however, we do not want to choose 9q such that the limited moment vector 9qU + (1 — 9q)uq lies 
exactly on the boundary of the realizable set 972.®, but rather so that the limited moment vector is in the 


Here we are using 1/2 instead of 1 to ensure that u and Ug are in the interior of 'R-Sl 

\uo<l 


13 





interior of 7^®. Therefore we define a tolerance e to add to each relevant Oqi (that is in [0,1]) as well as those 
facets such that 9qi G [—e, 0], indicating that while is on the correct side of the half space, it is closer 
than £ to the facet. Keeping in mind that 9g should not exceed one, this gives 


0 , := 


0 

min{l, £ + max{6»,i : 0qi G [-£, 1]}} 


if there is no 9gi G [—£, 1], 
else. 


Finally, in the implemented version, for the cell we set 0 = max{0g}. Figure 2a shows an example for the full 
moment model with N = 2. 



(a) Limiter example for N = 2 full moments p-basis with 
riQ = 7 and, for simplicity of exposition, Uqo = mq = 1. 



(b) The number of facets of for a few M^v (for N £ 
{4,6,8}) and MMjv (for N G {2,3,4}) models. The 
number of facets of MM 2 and MM 3 are almost exactly 
the same. 


Figure 2 


The main drawback to this implementation of the limiter is that, as illustrated in Figure 2b the number of 
facets d grows rapidly both with the number of moments N and the number of quadrature points Q. The 
numbers for this figure were computed using results from the study of convex polytopes, and a more detailed 
discussion is in the appendix. This issue was not a significant obstacle for our simulations but will be in 
higher dimensions where a much larger number of quadrature points is necessary. One clear speed up is 
that not every 9gi should be computed, as many facets do not intersect the line between u and u,j. More 
importantly, however, an implementation in higher-dimensions will probably have to further approximate 
7?.® by removing some irrelevant facets. 

Remark 3.4. There is actually no loss in accuracy by limiting to 7^® instead of when the quadrature Q 
defining 7?.® is the same as the quadrature used to compute all angular integrals: In this case we are solving 
the moment equations (2.41 with quadrature approximations to the angular integrals, and the numerical 
solution indeed lives on 7^®. (Furthermore, the function f can only be evaluated on 72.®.J The numerical 
solutions using quadrature then converge to moment equations with exact integrals as the quadrature rule 
converges. 


3.5. Slope limiting 

Although the realizability limiter plays some role in dampening spurious oscillations in numerical solutions, 
further dampening is needed [23j . We use the standard TVBM corrected minmod limiter proposed in [0]. 
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Assuming that the major part of the spurious oscillations are generated in the linear part of the underlying 
polynomial, whose slope in the j-th cell is simply uj, a basic limiter can be defined as 


^scalar ( 


u) = 


(^°) 

( 0 , 0 ,..., 0 ) 


. |uj| > M{/S.XjY and, 
otherwise 


for the j-th cell and the case fc = 2, that is piecewise quadratic approximations, so that the hnal row of 
zeros in the first case indicates that the coefficients for the quadratic basis functions are set to zero for each 
moment component. The absolute value and inequality are applied component-wise, and 


Dj{u) = D(Uj_i,Uj,Uj+i) := m(u],u°_^i - u°,u° - u°_i). 


The label “scalar” is used because the limiter is directly applied to each scalar component of u/j. The function 
m is the standard minmod function applied component-wise: 


TO(ai, 02, as) 


sign(ai)min{|ai|, |a 2 |, las]} if sign(ai) = sign(a 2 ) = sign(a 3 ), 
0 else. 


The constant M is a problem-dependent estimate of the second derivative, though we note that in [6] the 
authors did not find the solutions very sensitive to the value chosen for this parameter. 


However, it has been found that applying the limiter to the components themselves may introduce non¬ 
physical oscillations around an otherwise monotonic solution [5]. Therefore we instead apply the limiter to 
the local characteristic fields of the solution. The characteristic fields are found by transforming the moment 
vector u using the matrix V^, whose columns hold the eigenvectors of the Jacobian df/du evaluated at the 
cell mean Uj. We then transform back to the moment variables after applying the limiter. In the end, since 
Uj is a matrix of size (A: -|- 1) x {N + 1), this transformation is accomplished by post-multiplying with 
so that 

A,(u) = Af"i"’'(uV-^)Vj, 

where uV“^ is understood as (..., Uj_iV“^, UjV“^, Uj_|_iV“^,...). We apply this limiter to every Runge- 
Kutta stage. 


The Jacobian is computed at the cell means Uj using (3.4). This indeed also implies that we must solve the 
dual problem (2.6) to compute Q;(uj) for each cell. For cases when the matrix Vj has a condition number 
(defined using the two-norm) greater than a tolerance Atjac, we apply the scalar limiter. 


There are more advanced limiter strategies, for example WENO limiting |25l I38j or generalized minmod- 
limiting m, removing the drawback of having a problem-dependent parameter M. The slope limiter is 
however not a focus of this work. 


4. Numerical results 

We used the following parameter values: 
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r = 10 ® 

{r„} = {0,10-8,10-6,10-4,10-3} 

= 50 


uq = 40 
M = 50 


£= 10-44 

Kjac = lO® 


Optimization gradient tolerance. 

Outer regularization loop in optimizer. 

Number of optimization iterations before 
advancing outer regularization loop, 

Number of angular quadrature points. 

Slope-limiter constant, value suggested 
in [B], 

Realizability limiter tolerance. 

Condition-number tolerance for applying the characteristic 
transformation in the slope limiter. 


For the angular quadrature we used (nQ/2)-point Gauss-Lobatto rules over both ^ S [—1)0] and ^ G [0,1]. 
At the first time step, the initial multipliers for the optimizer (recall Section 3.3) are those of the isotropic 
density with the appropriate zeroth-order moment. For the rest of the simulation, the initial multipliers are 
those from the same point in space at the previous time step. 


To set the time step we use condition (3.13) with equality. 


4-1. Convergence tests 

We numerically test the convergence of the scheme two ways: First, we use the method of manufactured 
solutions on the total scheme, and second we focus on the convergence of the spatial reconstructions using 
the realizability limiter near the boundary of realizability. 


4.1.1. Manufactured solution 

In general, analytical solutions for minimum-entropy models are not known. Therefore, to test the convergence 
and efficiency of our scheme, we use the method of manufactured solutions. To avoid the effects of the 
boundary we use periodic boundary conditions, and we set the spatial domain to A = (—7r,7r). 

We begin by defining a kinetic density in the form of the entropy ansatz and which is periodic in space for 
every t: 


(j){t,x,fi) =ex.p{ao{t,x) + ai{t,x)fj,), 
a^ft, x) = — K — sin(a: — t) + c^t — ci, 
ai{t, x) =K + sin(a: — t). 

A source term is defined by applying the transport operator to (f: 

S(t, X, fj,) := dt4>(t, X, fj.) + gidx(j){t, x, fa). 


(4.1a) 

(4.1b) 

(4.1c) 


Thus by inserting this S into (2.1) (and taking Ua = CTs = 0) we have that (j) is, by construction, a solution of 

( 13 . 

A tedious but straightforward computation shows that choosing cq = 4 gives S > 0, which means that 
Theorem 3.3 will apply to the resulting moment system (for any K). Furthermore we take 


SO that the maximum value of {(f) for (t, x) G [0, t{] x A is one. The parameter K can be increased to make 
(j) look increasingly like a Dirac delta at ^ = 1. 
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Since our solution has the form of an entropy ansatz, v = (b0) is also a solution of (2.4) whenever 1 and ^ 
are in the linear span of the basis functions b. Clearly this holds for the M^r and MMat models for iV > 1. 
Notice also that v approaches the boundary of realizability as K is increased. 


We used the final time t{ = tt/S and chose K = 55, for which the maximum value of uijuQ is about 0.98 
(recall that |ui/uo| < 1 is necessary for realizability). In the following, we used the M 3 model so that our 
results included the effects of the numerical optimization. 

We compute errors in the zero-th moment of the solution, which we denote VQ(t,x) = Then 

and L°° errors for the zero-th moment uo^h{t,x) (that is, the zero-th component of a numerical solution Uh) 
are defined as 


E 


1 _ 

h 



uo,h{U,x)\dx and 


ma.x\vo{t{,x) - uo,h{tux )\, 

xGX 


(4.2) 


respectively. We approximate the integral in Ej^ using a 100-point Gauss-Lobatto quadrature rule over each 
spatial cell Ij, and E'^ is approximated by taking the maximum over these quadrature nodes. The observed 
convergence order v is defined by 


K2 \^X2j 


(4.3) 


where for i € { 1 , 2 }, Ej^^ is the error Ej^ for the numerical solution using cell size Aa;^, for p G { 1 , 00 }. 

A convergence table is presented in Table [fusing a tight gradient tolerance in the optimization of r = 10“^^. 
We observe that the expected convergence rates are achieved both in L^- and L°°-errors, although for fc = 0, 
the solution has only just begun to reach the convergent regime. In Figure [3a] we plot the L°°-error versus 
the computation time for the solution for the same value of r. Here we clearly see that higher-order methods 
are more efficient. 


The scheme is not convergent for arbitrarily large values of K. For large values of AT, the numerical solution 
will veer so close to the boundary of the realizable set that the optimization will have to use regularization, 
thus introducing errors into the solution. This was observed in [T], though here we can display this effect 
more precisely. In Figure [3bl we show the results using K = IIO for three spatial discretizations. In the 
most coarse discretization (J = 40), regularization is never necessary, bnt after doubling the number of cells, 
a few regularizations are used and their effects can be seen in the figure. Here, the optimizer regularizes 
four problems with r = 10“®, around x = —1.76, —1.68, 0.51, 0.74 at t = 0.38, 0.44, 0.35, 0.52, respectively, 
and the effect on the error spreads, mostly (to the right, the propagation direction of the solution), and 
magnifies slightly until the final time. The observed convergence order for Ej^ from J = 40 to 80 is = 2.9 
while the corresponding E^ order is = 2.1. When doubling the number of cells again (up to J = 160), the 
optimizer must regularize now only three problems with r = 10“®, this time around x = —1.74, —1.66, 1.04 
at t = 0.56, 0.51, 0.44, respectively. Indeed, immediately from the figure one can see that convergence in the 
E^ error has stopped and the observed convergence rate in E\ from J = 80 to 160 is only v = 1.6. 


4.I.2. Convergence of the limited spatial reconstructions 


Despite choosing a manufactured solution which lies close to the boundary of realizability, the realizability 
limiter was never active in any of our simulations in Section [4.1.11 Therefore in this section we artificially define 
a curve of moment vectors in space, reconstruct this curve in the finite-element space of discontinuous 
quadratic polynomials (recall (3.2)), and apply the limiter to move the reconstruction back into the set of 
numerically realizable moments. We then measure the convergence of this limited reconstruction!^ 


For simplicity, we only consider the monomial basis b = p. Using the Dirac delta function <5 = S{p,), we 
first choose two moment vectors Uq and Ui which lie arbitrarily close to the boundary of the numerically 


® We would like to thank an anonymous reviewer for suggesting such a test. 
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J 

o ■ 

II , 


k = 1 


k = 2 


K 

V 

El 

V 

El 

V 

20 

4.087e-h00 

— 

3.524e-02 

— 

1.897e-05 

— 

40 

2.608e-h00 

0.6 

9.532e-03 

1.9 

2.416e-06 

3.0 

80 

1.507e-h00 

0.8 

2.482e-03 

1.9 

3.049e-07 

3.0 

160 

8.161e-01 

0.9 

6.333e-04 

2.0 

3.828e-08 

3.0 

320 

4.256e-01 

0.9 

1.600e-04 

2.0 

4.796e-09 

3.0 

640 

2.174e-01 

1.0 

4.020e-05 

2.0 

6.072e-10 

3.0 


77'OO 

V 

77'OO 

Eh 

V 

77'OO 

Eh 

V 

20 

3.339e-01 

— 

3.164e-03 

— 

9.283e-06 

— 

40 

2.129e-01 

0.6 

8.543e-04 

1.9 

1.241e-06 

2.9 

80 

1.230e-01 

0.8 

2.222e-04 

1.9 

1.609e-07 

2.9 

160 

6.662e-02 

0.9 

5.666e-05 

2.0 

2.048e-08 

3.0 

320 

3.474e-02 

0.9 

1.431e-05 

2.0 

2.589e-09 

3.0 

640 

1.775e-02 

1.0 

3.595e-06 

2.0 

3.365e-10 

2.9 


Table 1: L^- and L°° -errors and observed convergence order u for the manufactured solution | |4.1| l with optimization gradient 
tolerance r = 10“^^. 


k — Q k — 1 k — 2 
leO 

le-2 

le-4 

le—6 

le-8 

le-10 

leO lei le2 le3 
Computation time (s) 

(a) L°°-Efficiency, r = 10“^^. Here we used the same 
numbers of cells as in Table 


^ J = 40 J = 80 ^ J = 160 



X 

(b) Regularization destroying convergence for K = 110. 



Figure 3: Using the manufactured solution to consider the efficiency of higher-order methods and the effects of regularization. 
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Uo Ui 


J 

K 


77 ’00 

u 

K 

u 

77 ’00 

ly 

^max 

8 

1.072e-03 

— 

1.848e-03 

— 

5.790e-04 

— 

9.983e-04 

— 

1.56e-02 

16 

1.130e-04 

3.2 

2.469e-04 

2.9 

6.108e-05 

3.2 

1.334e-04 

2.9 

4.27e-03 

32 

1.342e-05 

3.1 

3.137e-05 

3.0 

7.253e-06 

3.1 

1.695e-05 

3.0 

1.39e-03 

64 

1.655e-06 

3.0 

3.937e-06 

3.0 

8.943e-07 

3.0 

2.127e-06 

3.0 

6.61e-04 

128 

2.060e-07 

3.0 

4.927e-07 

3.0 

1.113e-07 

3.0 

2.662e-07 

3.0 

4.48e-04 

256 

2.564e-08 

3.0 

6.160e-08 

3.0 

1.385e-08 

3.0 

3.328e-08 

3.0 

2.51e-04 

512 

3.188e-09 

3.0 

7.700e-09 

3.0 

1.723e-09 

3.0 

4.160e-09 

3.0 

3.76e-06 


Table 2: L^- and L°°-errors and observed convergence order u for the first two components of the realizability-limited, piecewise 
quadratic reconstruction of u(a:) from \4A\ with 7 = 10“^*^ and moment order N = 4. 


realizable set: 


Uo := (1 - 7 ) {p6{fi - Ho)) + 7Uiso = (1 - 7 )p(Mo) + 7Uiso 

ui := 10“® ((1 - 7 ) {pS{p, + 1)) + 7 Uiso) = 10“® ((1 - 7)P(-1) + 7 Uiso) 


where ho ~ 0.5403, which is a node from the angular quadrature we use (see the beginning of Section |^, Ujgo 
are the moments of the isotropic density (see (3.9|), and 7 G [0,1] controls the distance to the boundary. For 
iV > 1, both Uq and Ui lie on the boundary of the realizable set when 7 = 0 . Since the Dirac delta functions 
are placed at nodes in the angular quadrature, Uq and Ui (and any convex combination thereof) are in 7^® 
for 7 G (0,1], and so we define a curve of moments in space by taking convex combinations of Uq and Ui: 


u(a:) := (1 - A(a;))uo + A(a;)ui, x G [-1,1], 


(4.4) 


where A(a;) € [0,1] is given by 


A(q 


cos(7ra;) + 1 
2 


X G [-1,1]. 


To perform the convergence test, we project u(a:) onto for increasing numbers of cells J and then apply the 
realizability limiter. Errors and observed convergence order are computed as in (4.21 and (4.3), respectively, 
over the finest grid of the tests, here J = 512 cells. 


We found that taking 7 G [10“^°, 10“^] places the moment curve u(a:) close enough to the boundary of 
realizability that the realizability limiter was active for every number of cells we considered. In Table 
we show convergence rates for the smallest 7 in this range. These results show the desired third-order 
convergence. In this table we include the column ^max, which gives the maximum value of 9 from the 
realizability limiter over all spatial cells. That dmax is nonzero in each row indicates that the realizability 
limiter was active for every reconstruction. We observed similar results for every moment component and 
moment curves of any order N G {2, 3,4, 5}. 


However, it has been remarked in [36] that in some pathological situations this limiter may reduce accuracy 
to second order. We can demonstrate this by pushing u(a;) closer to the boundary of realizability by setting 
7 = 10“^^. These results are given in Table where in the error we still see third-order convergence, but 
in the error the convergence has degraded to second order. 


In further tests, which for brevity we omit here, we observed similar behavior as we decrease 7 until we 
arrive at 10“^"*^, which is the value of the tolerance in the realizability limiter (see the beginning of Section]^. 
At this point convergence is only first order as expected. 


That the limiter works so close to the boundary of realizability is satisfactory to us: When moment vectors as 
close as 10 “^° to the boundary of realizability appear in a simulation, errors from the numerical optimization 
are likely to affect convergence before these effects of the realizability limiter play a role. 
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Uo Ui 


J 

K 


77’00 

u 

K 

u 

77’00 

ly 

^max 

8 

1.308e-03 

— 

1.848e-03 

— 

7.066e-04 

— 

9.983e-04 

— 

1.93e-02 

16 

1.224e-04 

3.4 

2.469e-04 

2.9 

6.614e-05 

3.4 

1.334e-04 

2.9 

8.11e-03 

32 

1.460e-05 

3.1 

3.137e-05 

3.0 

7.889e-06 

3.1 

1.695e-05 

3.0 

5.23e-03 

64 

1.805e-06 

3.0 

7.105e-06 

2.1 

9.754e-07 

3.0 

3.839e-06 

2.1 

4.55e-03 

128 

2.247e-07 

3.0 

1.721e-06 

2.0 

1.214e-07 

3.0 

9.299e-07 

2.0 

4.32e-03 

256 

2.799e-08 

3.0 

4.162e-07 

2.0 

1.512e-08 

3.0 

2.248e-07 

2.0 

4.15e-03 

512 

3.457e-09 

3.0 

8.958e-08 

2.2 

1.868e-09 

3.0 

4.840e-08 

2.2 

3.57e-03 


Table 3: L^- and L°°-errors and observed convergence order u for the first two components of the realizability-limited, piecewise 
quadratic reconstruction of u(a:) from \AA\ with 7 = 10“^^ and moment order N = i. 


4-2. Plane source 

In this test case we start with an isotropic distribution where the initial mass is concentrated in the middle 
of an infinite domain x S (— 00 , 00 ): 


tpt=oix, /t) = V'floor + S{x) 

where the small parameter ipdoor = 0.5 x 10“® is used to model a vacuumj^ In practice, a bounded domain 
must be used, so we choose a domain large enough that the boundary should have only negligible effects on 
the solution: thus for our final time tf = 1, we take X = [a:L,iCR] = [—1-2,1.2]. At the boundary we set 

n) = IpHoor and M) = V'floor 


We use isotropic scattering and no absorption, therefore C {ijj) = ^{tp) — tp, as = 1 and Ua = 0. 

We approximate the delta function by using an even number of spatial cells and splitting the delta into the 
cells immediately to the left and right of a; = 0. This is then projected into V/p using (3.3b) ^ All solutions 
here are computed with J = 300 cells and spatial polynomials of degree k = 2. 


Figure [^presents solutions for different moment models. This figure includes a reference solution, which we 
computed using our scheme for the Pgg model with J = 2000 cells and spatial order k = 0. The oscillations 
we see have been observed before and arise due to the fact that we are using moment models (which are 
indeed spectral methods) on a non-smooth problem. Indeed, our M^v results agree well with those in |I2j . 
The full- and mixed-moment models look largely similar for the same numbers of degrees of freedom, with 
the notable differences being around x = 0. Here the mixed-moment solutions are much more sharply peaked 
while the full-moment solutions are more flat and wider. The discrepancy in magnitude, however, seems to 
decrease as N increases, which agrees with the expectation that both methods are converging as N ^ 00 . 


A vacuum is not exactly realizable by the entropy ansatz 113 for the Maxwell-Boltzmann entropy. 

We refer the reader to m for a thorough discussion of using discontinuous Galerkin methods for problems with delta 
functions in the initial data. 
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Figure 4: Local density uq for different models and a reference solution at t = 0.8 in the plane-source problem. 


Figure shows the activity of the realizability limiter. The realizability limiter is most active along the 
fronts where particles from the initial impulse are first entering the domain, which is as expected, because 
this is where the moment vectors of the solution lie closest to the boundary of realizability [2] . We also see 
that as the number of moments increases, the activity of the realizability limiter increases as well. This is 
also not surprising, since realizability conditions typically require tighter and tighter bounds on the moment 
components as their order increases. Thus numerical errors of a similar size will have an increasingly large 
chance of pushing the solution out of the realizable set. We did not observe significant differences between 
the full-moment and mixed-moment models. 
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Figure 5: The value of 9 in the realizability limiter for two models of the plane-source problem. Note that we choose a 
logarithmic scale so that even small values of 9 are noticeable. 


4-3. Two beams 

This test case models two beams entering a (nearly) empty absorbing medium. This is a classical test problem 
used to illustrate the shortcomings of the Mi model, whose steady-state solution for this problem has a 
nonphysical shock. This test case is also challenging for the numerical optimization [5]. 
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The domain is X = (—0.5,0.5), and the beams are specified by boundary conditions 


= ^exp ( - 


2E2 


and fJ.) = ^ exp [- 


2E2 


with E = 50. The initial condition again approximates a vacuum similarly as before: 

tjj{0,x,p,) = '(/jfloor = 0.5 X 10“®. 


Finally, we use absorption parameter CTq = 2 and no scattering, CTs = 0. Again, all solutions here are computed 
with J = 300 cells and spatial polynomials of degree k = 2. 

Our numerical solutions for the local density Uq all look qualitatively the same—with the notable exception 
of the Ml model—and match the true steady-state solution, so in Figure [6a| we only present one example 
solution for the unfamiliar reader. We do note, however, that our steady-state solutions for the full-moment 
models do not contain the shocks observed in [12]. This difference is apparently due to the fact that we use 
sharply forward-peaked boundary conditions, where as isotropic boundary conditions were used in |12j . As 
shown in Figure [ba] the mixed-moment solutions do not appear to contain any steady-state shocks either. 

In Figure [6b| we present a zoomed-in view of the oscillations present in a typical transient solution. These 
oscillations become more noticeable for higher-order models (M^r for > 4 and MM^r for N > 3). It 
seems clear to us that these are numerical artifacts, and we believe they would be mitigated with a more 
sophisticated slope limiter. 




(b) The M 4 model at t = 0.8. 


Figure 6 : Local density uq for different models at t = 0.8 in the two-beam problem. 


The activity of the realizability limiter again increases with the number of moments, but in this problem we 
see differences between full- and mixed-moment models. Figure [Tj illustrates this difference. The reason for 
this difference is not yet clear to us, but it seems to indicate that the mixed-moment model is converging 
more slowly to steady state. 


5. Conclusions and outlook 

We presented a high-order Runge-Kutta discontinuous Galerkin scheme for minimum-entropy moment models 
of linear kinetic equations in one space dimension. The key issue for higher-order methods for minimum- 
entropy moment models is that the numerical solution typically leaves the set of realizable moments, even 
though standard techniques can be used to show that the cell means of the solution remain realizable. We 
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Figure 7: The value of 9 in the realizability limiter for two models of the two-beam problem. Note that we choose a logarithmic 
scale so that even small values of 9 are noticeable. 


address this problem using a realizability limiter inspired by the positivity-preserving limiter used in |35j 
for the Euler equations. Such a limiter requires the computation of the intersection of a line in moment 
space with the boundary of the realizable set, a set which typically has nonlinear boundaries. We are able to 
approximate this intersection by replacing the true realizable set with its quadrature-based approximation, 
which is a convex polytope. This quadrature-based approximation is intriguing because it is a convex polytope 
for any moment order and any dimension of the angular domain indicating that our techniques could be 
extended to these cases. 

We constructed a new manufactured solution whose source term is realizable, thus allowing us to consider 
target solutions closer to the boundary of the realizable set. These tests show that our scheme converges as 
expected and that higher-order schemes are more efficient. We also present numerical solutions for standard 
benchmark problems, where we are able to compare full- and mixed-moment models. 

Future work should focus on a parallelized implementation for two- and three-dimensional problems. Theo¬ 
retically, implementation of the quadrature-based realizability limiter requires no change because the convex 
polytopic structure of 7^® holds in any dimension. Practically speaking, however the main challenge is that 
the number of facets grows quickly with the number of moments and number of quadrature points, both of 
which will be higher. Further work in higher-order methods will also have to consider new methods for time 
integration, as here we relied heavily on the SSP property, which is not possible past fourth-order. Relatedly, 
at least partially implicit time integrators should be investigated, particularly in the context of constructing 
an asymptotic preserving numerical method for the moment system. 


Appendix A. The number of facets in 7?.® 

Even with some speed-ups in the computation of the facet-intersections and possibly approximations by 
removing facets, the number of facets plays a large role in determining the complexity of finding the 
intersection of a ray with the boundary of the convex polytope In this section we mention how 

some results from the study of convex polytopes give the exact number of facets in the full-moment case, 
and then we compare this with an upper bound of the number of facets in the mixed-moment case. 

First, some notational remarks for this section: For convenience, we work with the closures 72.®|„g<i and 
of Lo<i Lp-i respectively. When working with 7?.®|„g=i, we consider it as a subset of 


and n^\ 
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(or in the mixed-moment case), and use the notation Ui and bi to indicate the final N (or 2N) 
entries of u and b respectively. We also often work with the matrix form of the half-space representation, 
so for example = {ui : ylui < 6}, for a matrix A G with rows and a vector b G 

Finally, we omit proofs in this section because we consider the arguments needed to be unenlightening and 
relatively straightforward. 


We first note that the number of facets of 7^®l“o<i more than that of 

Proposition Appendix A.l. If A and b define a half-space representation of TZ^\u„=i, then a half-space 
representation of Tl^\ug<i is: 



Uo<l 



( Uq 

Vui 




Therefore in the sequel we focus on the number of facets of 72.®|„o=i. 

First we consider the full-moment case. The convex polytope TZ^\ua=i C is known as the cyclic polytope 
and plays a special role in the study of convex polytopes. The Upper Bound Theorem states that for a given 
number of vertices in a given dimension, the cyclic polytope has the maximum number of facets [1]. Gale’s 
evenness condition or the Dehn-Sommerville equations can be used to show that the number of facets is 

C(JV,„e) = ("0 - [|(iV+ 1)J\ ^ fna - [1(JV + 2)J 
V nQ-N J \ uq-N 

for UQ > A' > 1 dj, where [-J indicates the integer part of its argument. We note that this holds for any 
choice of distinct quadrature nodes {^i}. Since 72.p|„o=i has C{N,nQ) facets, there exists a half-space 
representation such that A G and b G Unpacking the definition of the binomial 

coefficient we can see that for fixed, even N, we have C{N,nQ) = and for fixed, odd N we have 

C(A,ne) = 0(ne(^-i)/2). 

One can see from Figure !^ that MMjv models appear always to have fewer facets than the corresponding 
Mjv models. To show that this holds more generally, we first need a half-space representation for 7?.mUo=i- 
This representation can be derived using the half-space representations from the full-moment case. 

Proposition Appendix A.2. Let A± and b± define half-space representations for the convex polytopes 
formed by the basis functions on the positive and negative subintervals respectively: 


co{pi(/i*)}^^>o = {ui+ : A+U1+ < &+}, 
co{pi(m*)}^_<o = {ui- : ^-ui- < 

We assume b± > 0 component-wise^^ Then a half-space representation for TZm\uo=i given by 

= (nf-) ■ ^ ‘I ■ 

where 



/ 



fbA 

A = 

0 

A_ 

and b = 

b_ 


[ ■ 

■■ J 


1 

1:/ 


^^Here we are using the fact that the subintervals for the mixed-moments are joined exactly at ^ = 0 and assume furthermore 
that this point is a quadrature node. This is indeed a reasonable assumption, since even in MMi, a delta function can form at 
fj, = 0. 
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where the last rows of the A include only those pairs {i,j} such that neither b^i nor b-j is equal to zero. 

If we let C± denote the number of rows of A± respectively, this representation gives C+ + C- + C+C- as 

an upper-bound on the number of facets in 72.mUo=ilH] The number of rows such that b±i = 0 is equal to 
the number of facets including the vertex corresponding to the quadrature point at /i = 0. These facets 
can be more generally described as those containing the vertex corresponding to the first quadrature point, 
when the quadrature points are arranged in increasing order. The number of such facets can be computed 
using Gale’s evenness condition (see Theorem 13.6 and Exercise 13.1 in [3]). We omit this computation here 
but note that removing these facets does not change the order of the number of facets (nor any relevant 
leading-order coefficients), so in the comparison that follows, we ignore these terms. 

To compare the full-moment and mixed-moment cases for the same number of degrees of freedom, one would 
consider the full-moment case of order iV, for N even, and the mixed-moment case of order N/2. Let us 
assume that we use a quadrature set which includes fj, — 0 and has Q/2 points over both /i > 0 and fJ. < 0, 
for a total of Q — 1 points (since the point at /i = 0 should not be counted twice in the full-moment case). 
Then the number of facets in the full-moment case is C{N,Q — 1) while in the mixed-moment case, the 
number of facets in our half-space representation is on the order of C'(iV/2, <5/2)^. Then, straightforward 
calculations show that when N/2 is odd we have C{N/2, Q/2)'^ = which is one order less than 

in the full-moment case. When N/2 is even, our half-space representation for the mixed-moment case has 
facets, which is the same order as the full-moment case. However, the leading-order coefficient is 
smaller in the mixed-moment case, thereby showing that the number of facets in the mixed-moment case 
is at least asymptotically smaller. Indeed, if we let N = 4n, the ratio of the highest-order coefficients is 
(2n)!/(2”n!)^, which is bounded by 1/2 and monotonically decreases with n. 
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